﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;
using dnAnalytics.Math;
using CalcSharp.Core.Properties;

namespace dnAnalytics.Statistics.Distributions
{
    /// <summary>
    /// The Gamma distribution is a distribution over the positive real numbers. It is parameterized by a shape
    /// parameter a > 0 and an inverse scale parameter b > 0. The pdf is
    /// 
    ///     p(x) = \frac{b^a}{\Gamma(a)} x^{a-1} \exp{-bx}.
    /// </summary>
    /// <remarks>The distribution will use the <see cref="System.Random"/> by default. 
    /// Users can set the random number generator by using the <see cref="RandomNumberGenerator"/> property.
    /// 
    /// The statistics classes will check all the incoming parameters whether they are in the allowed
    /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
    /// to false, all parameter checks can be turned off.</remarks>
    public class Gamma : IContinuousDistribution
    {
        // Gamma shape parameter.
        private readonly double mA;

        // Gamma scale parameter.
        private readonly double mB;

        /// <summary>
        /// Constructs a Gamma distribution.
        /// </summary>
        /// <param name="a">The shape of the Gamma distribution.</param>
        /// <param name="b">The inverse scale of the Gamma distribution.</param>
        public Gamma(double a, double b)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(a, b);
            }

            mA = a;
            mB = b;
            RandomNumberGenerator = new System.Random();
        }

        /// <summary>
        /// The Gamma shape parameter.
        /// </summary>
        public double A
        {
            get { return mA; }
        }

        /// <summary>
        /// The Gamma scale parameter.
        /// </summary>
        public double B
        {
            get { return mB; }
        }

        #region IDistribution Members

        /// <summary>
        /// The mean of the distribution.
        /// </summary>
        /// <value></value>
        public double Mean
        {
            get { return mA/mB; }
        }

        /// <summary>
        /// The standard deviation of the distribution.
        /// </summary>
        /// <value></value>
        public double StdDev
        {
            get { return System.Math.Sqrt(mA/(mB*mB)); }
        }

        /// <summary>
        /// The variance of the distribution.
        /// </summary>
        /// <value></value>
        public double Variance
        {
            get { return mA/(mB*mB); }
        }

        /// <summary>
        /// The entropy of the distribution.
        /// </summary>
        /// <value></value>
        public double Entropy
        {
            get { return mA - System.Math.Log(mB) + SpecialFunctions.GammaLn(mA) + (1.0 - mA)*SpecialFunctions.DiGamma(mA); }
        }

        /// <summary>
        /// Gets or sets the random number generator.
        /// </summary>
        /// <value>The random number generator used to generate a random sample.</value>
        public System.Random RandomNumberGenerator { get; set; }

        #endregion

        #region IContinuousSamplingDistribution Members

        /// <summary>
        /// The mode of the distribution.
        /// </summary>
        /// <value></value>
        public double Mode
        {
            get { return (mA - 1.0)/mB; }
        }

        /// <summary>
        /// The median of the distribution.
        /// </summary>
        /// <value></value>
        public double Median
        {
            get { throw new Exception("Not implemented yet."); }
        }

        /// <summary>
        /// Evaluates the probability density function for the Gamma distribution.
        /// </summary>
        /// <param name="x">The location in the domain where we want to evaluate the probability density function.</param>
        /// <returns></returns>
        public double Density(double x)
        {
            return System.Math.Pow(mB, mA)*System.Math.Pow(x, mA - 1.0)*System.Math.Exp(-mB*x)/SpecialFunctions.Gamma(mA);
        }

        /// <summary>
        /// Samples a Gamma distributed random variable.
        /// </summary>
        public double Sample()
        {
            return SampleGamma(RandomNumberGenerator, mA, mB);
        }

        /// <summary>
        /// Samples an array of Gamma distributed random variables.
        /// </summary>
        /// <param name="n">The number of variables needed.</param>
        public double[] Sample(int n)
        {
            return DoSample(RandomNumberGenerator, n, mA, mB);
        }

        #endregion

        /// <summary>
        /// Evaluates the cumulative distribution function for the Gamma distribution.
        /// </summary>
        public double CumulativeDistribution(double x)
        {
            return SpecialFunctions.IncompleteGamma(mA, x*mB, true);
        }

        /// <summary>
        /// Evaluates the logarithm of the probability density function for the Gamma distribution.
        /// </summary>
        public double DensityLn(double x)
        {
            return mA*System.Math.Log(mB) + (mA - 1.0)*System.Math.Log(x) - mB*x - SpecialFunctions.GammaLn(mA);
        }

        /// <summary>
        /// Sampling implementation based on:
        /// "A Simple Method for Generating Gamma Variables" - Marsaglia &amp; Tsang
        /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="alpha">The shape of the Gamma distribution.</param>
        /// <param name="beta">The inverse scale of the Gamma distribution.</param>
        /// <returns></returns>
        protected static double SampleGamma(System.Random rnd, double alpha, double beta)
        {
            double a = alpha;
            double alphafix = 1.0;

            // Fix when alpha is less than one.
            if (alpha < 1.0)
            {
                a = alpha + 1.0;
                alphafix = System.Math.Pow(rnd.NextDouble(), 1.0/alpha);
            }

            double d = a - 1.0/3.0;
            double c = 1.0/System.Math.Sqrt(9.0*d);
            while (true)
            {
                double x = Normal.Sample(rnd, 0.0, 1.0);
                double v = 1.0 + c*x;
                while (v <= 0.0)
                {
                    x = Normal.Sample(rnd, 0.0, 1.0);
                    v = 1.0 + c*x;
                }
                v = v*v*v;
                double u = rnd.NextDouble();
                x = x*x;
                if (u < 1.0 - 0.0331*x*x)
                {
                    return alphafix*d*v/beta;
                }
                if (System.Math.Log(u) < 0.5*x + d*(1.0 - v + System.Math.Log(v)))
                {
                    return alphafix*d*v/beta;
                }
            }
        }

        /// <summary>
        /// Samples a Gamma distributed random variable.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="a">The shape of the Gamma distribution.</param>
        /// <param name="b">The inverse scale of the Gamma distribution.</param>
        public static double Sample(System.Random rnd, double a, double b)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(a, b);
            }

            return SampleGamma(rnd, a, b);
        }

        /// <summary>
        /// Samples an array of Gamma distributed random variables.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of variables needed.</param>
        /// <param name="a">The shape of the Gamma distribution.</param>
        /// <param name="b">The inverse scale of the Gamma distribution.</param>
        public static double[] Sample(System.Random rnd, int n, double a, double b)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(a, b);
            }

            return DoSample(rnd, n, a, b);
        }

        /// <summary>
        /// Check the parameters of the Gamma distribution.
        /// </summary>
        /// <param name="a">The shape of the Gamma distribution.</param>
        /// <param name="b">The inverse scale of the Gamma distribution.</param>
        /// <exception cref="ArgumentOutOfRangeException">If any of the Gamma parameters are non-positive.</exception>
        private static void CheckParameters(double a, double b)
        {
            if (a <= 0.0)
            {
                throw new ArgumentException("a", Resources.NotPositive);
            }
            if (b <= 0.0)
            {
                throw new ArgumentException("b", Resources.NotPositive);
            }
        }


        /// <summary>
        /// Samples an array of Gamma distributed random variables.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of variables needed.</param>
        /// <param name="a">The shape of the Gamma distribution.</param>
        /// <param name="b">The inverse scale of the Gamma distribution.</param>
        private static double[] DoSample(System.Random rnd, int n, double a, double b)
        {
            double[] arr = new double[n];
            for (int i = 0; i < n; i++)
            {
                arr[i] = SampleGamma(rnd, a, b);
            }
            return arr;
        }
    }
}